On the spectrum of fluctuations of a liquid surface: 
Prom the molecular scale to the macroscopic scale 
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We show that to account for the full spectrum of surface fluctuations from low scattering vector 
qd<^ \ (classical capillary wave theory) to high qd>l (bulk- like fluctuations), one must take account 
of the interface's bending rigidity at intermediate scattering vector qd< 1, where d is the molecular 
diameter. A molecular model is presented to describe the bending correction to the capillary wave 
model for short-ranged and long-ranged interactions between molecules. We find that the bending 
rigidity is negative when the Gibbs equimolar surface is used to define the location of the fiuctuating 
interface and that on approach to the critical point it vanishes proportionally to the interfacial 
tension. Both features are in agreement with Monte Carlo simulations of a phase-separated colloid- 
polymer system. 



I. INTRODUCTION 

The description of the spectrum of surface 
fluctuations of a liquid from the macroscopic 
scale down to the molecular scale remains a chal- 
lenging experimental and theoretical problem. 
Using grazing incidence light scattering exper- 
iments, Daillant and coworkers^ were able, for 
the first time, to determine the full spectrum of 
surface fluctuations, where in previous experi- 
ments (ellipsometry, reflectivity) only certain as- 
pects of the spectrum could be determined. At 
the same time, the spectrum can now be ana- 
lyzed in computer simulations with ever increas- 
ing accuracy^i^i^. 

Theoretical insight into the structure of a sim- 
ple liquid surface is provided by density func- 
tional theories on the one hand^i^i and the 
capillary wave model on the other hand^iSiiiS. 
Density functional theories provide a descrip- 
tion of the interface on a microscopic level. 
The prototype of such theories, the van der 
Waals squared-gradient model, was very suc- 
cessful in describing, for the first term, the den- 
sity profile and surface tension in terms of molec- 
ular parameters^. It, however, fails to cap- 
ture the subtle role of long wavelength interfa- 
cial fiuctuations described by the capillary wave 
model^>a>ia. 

The capillary wave model introduced in 1965^ 
describes the spectrum of fiuctuations in terms 
of a height function /i(r||) with the surface ten- 
sion a and gravity g acting as the dominant 
restoring forces. The length scale involved in de- 
scribing capillary waves is the capillary length, 
Lc = a/ a/ (m Apg), which may be as large as a 
tenth of a millimeter. The theoretical challenge 
is to incorporate both theories and to describe 
the spectrum of fiuctuations of a liquid surface, 
as determined from light scattering experiments 
and computer simulations, from the molecular 
scale to the scale of capillary waves. 

An important ingredient in "bridging the 



gap" between capillary waves and the molecular 
scale is an extension of the capillary wave model 
that incorporates the energy associated with 
bending the interfac o^^d^d^ . Bending is impor- 
tant when the wavelength of the height fiuctu- 
ations is approximately y^k-gT/a, which is typ- 
ically of the order of a few times the molecular 
diameter, i.e. close to the scale where the molec- 
ular structure becomes important and the den- 
sity fiuctuations are more bulk-like. The natural 
question that arises is whether it is possible to 
describe the full spectrum of surface fiuctuations 
by the capillary wave model at long wavelengths 
and bulk-like fluctuations at the molecular scale. 
Is it then necessary to include the leading or- 
der correction to the capillary wave model from 
bending or are even higher order terms, relevant 
at even smaller length scales, required? 

This article addresses these questions in two 
partsfa condensed version has appeared in 
ref. \14 ). In the flrst part, we analyze the 
spectrum of fiuctuations recently obtained by 
Vink et al.^ in computer simulations of a phase- 
separated polymer-colloid syste m-'^^d^d'^d^ in 
which the interactions are strictly short-ranged. 
It is shown that the simulation data are very 
accurately described by the combination of the 
capillary wave model extended to include a 
bending correction, with the bending rigidity as 
an adjustable parameter, and bulk-like fiuctua- 
tions. 

In the second part, a molecular basis for the 
bending correction to the capillary wave model 
is offered and the results are compared with the 
simulations. The theoretical framework used for 
the comparison is mean-field density functional 
theory in which the interactions are described by 
a non-local, integral terro^i^iiii^. The advantage 
of this approach is that it features the full shape 
of the interaction potential enabling the analysis 
of different forms and ranges of the interaction 
potential. We consider both short-ranged inter- 
actions, and long-ranged interactions, that fall 
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of as C/(r) oc 1/ at large intcrmolecular sepa- 
rations. 

An important ingredient in our theoretical 
analysis is the modification of the density pro- 
file, described by pi{z), due to the local bend- 
ing of the interface^*'. The determination of 
Pi{z) requires one to formulate precisely the 
thermodynamic conditions used to vary the in- 
terfacial curvature. Several approaches for the 
determination of pi{z) have appeared in the 
literaturei^i^i2ii^. They differ in the form of 
the external field used to set the curvature to 
a specific value; in the equilibrium approaches 
the external field is uniform throughout the sys- 
tem, whereas in the approach by Parry and 
Boulter— it is infinitely sharp-peaked (T^xt 
5{z)) at the interface. In this article we suggest 
to add an external field acting in the interfacial 
region only with a peak- width of the order of the 
thickness of the interfacial region. The advan- 
tage of this approach is that the bulk regions are 
unaffected by the additional of the external field 
and the resulting pi{z) is a continuous function. 

Our paper is organized as follows: in Section 
2, the general form of the surface structure fac- 
tor to describe the spectrum of interfacial fluctu- 
ations is derived as the combination of the cap- 
illary wave model extended to include a bend- 
ing correction and bulk-like fluctuations. This 
form is then compared in Section 3 to the Monte 
Carlo (MC) simulation results by Vink et al.^ 
for the phase-separated polymer-colloid system. 
In Section 4, the mean-field density functional 
theory used to provide a molecular basis for the 
bending extension to the capillary wave model 
is presented. Explicit results are obtained for 
short-ranged interactions (Section 5), and long- 
ranged interactions (Section 6). We end with a 
discussion of results. 



II. THE FLUCTUATING LIQUID 
SURFACE 

In the classical capillary wave model (CW), 
the fluctuating interface is described by a two- 
dimensional surface height function /i(?^|), where 
r|| = {x^y) is the direction parallel to the 
surface^iSiiiS. The fluctuating density profile can 
then be written in terms of an "intrinsic density 
profile" shifted over a distance h{r]\): 



Pir) = Pq{z - /i(f||)) 



(1) 



where po{z) is the intrinsic density profile. Of- 
ten, fluctuations are assumed to be small so that 
an expansion in h can be made, neglecting terms 
of 0(/i2), 



An important consequence of the above lin- 
earization is that one may now identify the in- 
trinsic density profile as the average density pro- 
file, po{z) =< p{r) >, in view of the fact that 
< h{r\\) >— 0. It is convenient to locate the 
z — Q plane such that it coincides with the Gibbs 
equimolar surface^i^, i.e. 



dz[<p{r)> -pstcpiz)] 

o 

= Idzlpo (z) - pstcp (^) ] = , (3) 



where PstGp(z) = Pe ©(--z) + Pv Qiz) with 9(z) 
the Heaviside function and pe^v the bulk density 
in the liquid and vapor region, respectively. 

In the above model for p(r), the density cor- 
relations are essentially given by the correla- 
tions of /i(r[| ), which are described by the height- 
height correlation function: 



Shhir\\) = </).(fi_[|)ft.(f2j|)> 



(4) 



where =r2,|l ^ and r|[ =| fj| |. 

To determine the height-height correlation 
function, one should examine the change in free 
energy. Ail, associated with a fiuctuation of the 
interface. In the capillary wave model it is de- 
scribed by considering the change in free en- 
ergy associated with a distortion of the surface 
against gravity and surface area extension^: 



Ar2 = i Idrii 



m Ap g h{r 



It is convenient to express Afi 
of the Fourier Transform of h(f\\ 



in terms 
hiq) = 



[mApg + aq'^] h{q)h{-cf). 



2 J (27r)2 

(6) 

In the capillary wave model the height-height 
correlation function is determined by a full Sta- 
tistical Mechanical analysis^^^ in which the 
above expression for the change in free energy 
is interpreted as the so-called capillary wave 
Hamiltonian, Afl — T-lcw[h{f\\)]. In general, one 
has 



hh 



I h{r2, 



(7) 

where Z is the partition function associated with 
'Hcw[h], fee is Boltzmann's constant and T is the 
temperature. It can be shown that^ii^ 



Shh{q) 



p{r) = pn{z)- p'Q{z)h{r\\) + 



dr\\ e-*«-'-|i^^;,(r||) 



(8) 



(2) 



mApg + aq^ a [Lc + q^) 
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For simplicity, we ignore gravity effects in the 
following and set Lc = oo (g = 0). 



A. Extended capillary wave model 

In the derivation of the classical capillary 
wave model, one assumes an expansion in gra- 
dients of h{r\\), I V/i |<C 1. In the extended cap- 
illary wave model (ECW), one wishes to extend 
the expansion by including higher derivatives of 
h{r\\). To leading order one may then write the 
fluctuating density asi^^ 



Ah{fi 



(9) 

The function pi{z) is identified as the correction 
to the density profile due to the curvature of the 
interface, A/i(r^|)«— l/i?i — I/R2, with i?i and 
i?2 the (principal) radii of curvature. The pref- 
actor of —1/2 is chosen such that the notation is 
consistent with an analysis in which the curva- 
ture does not result from a fluctuation of the pla- 
nar interface, but is due to the fact that one con- 
siders a spherical liquid droplet {Ri = i?2 = i?) 
in (metastable) equilibrium with a bulk vapor 
phase^i^ii^. An expansion in the curvature of 
the density profile Ps(r) then gives 



Psir) = Pair) + 



R 



(10) 



which parallels the expansion in Eq.Q. 

The inclusion of curvature corrections in the 
extended capillary wave model leads to higher 
order terms in an expansion in , terms beyond 
aq^, in the expression for Afl in Eq.®. It is 
customary to capture these higher order terms 
by introducing a wave vector dependent surface 
tension a{qy^ 



dq 
(2^ 



a{q)q'h{q)h{^q), (11) 



which gives for the height-height correlation 
function 



Shhiq) = 



kBT 



(12) 



The precise form of (T{q) depends sensitively on 
the behavior of the interaction potential at large 
distancesi^. When the interaction potential is 
sufficiently short-ranged (SR), the expansion of 
a{q) in q^ is regular and the leading correction 
is of the form: 



a{q)=a + kq^ +0{q^). 



(SR) (13) 



The coefficient k is identified as the bending 
rigidity^i^^^. This is because the form for 



An in Eq.dni), with (j{q) given by Eq.lUS]), can 
also be derived from the Helfrich free energy 
expressionii, which reads for a fluctuating in- 
terface: 



An 



dru 



CT|V/i(f||)P + fc(A/i(r||)) 



(14) 

When the interaction potential is long-ranged 
(LR), specifically when it falls of as U{r)(xl/r^ 
at large intermolecular distances, which is the 
case for regular fluids due to London-dispersion 
forces, one finds that the leading correction to 
(T(g) picks up a logarithmic contributior>i^ : 

aiq)=a + k,q^ln{qh)+Oiq^), (LR) (15) 

with ks and £k parameters independent of q. 
The coefficient ks depends on the asymptotic 
behavior of U (r) but is otherwise a universal 
constant^^. The bending length ik depends, 
like the bending rigidity k, on the microscopic 
parameters of the model. In principal, all the 
parameters tr, k, kg, and £k can be expressed 
in terms of the density profiles po{z) and pi{z) 
by inserting the fluctuating density as given in 
Eq.® into a microscopic model for the free en- 
ergy and comparing the result with Eq. (|lip . 

It is important to realize that the extended 
capillary wave model assumes a curvature ex- 
pansion in Eq.Q which translates into an ex- 
pansion in q^ in Eg. ljlip that is valid only up 
to 0{q'^). Higher order terms are not system- 
atically included. The result is that one should 
fimit the expansion of a{q) in Ea.([T3|) or Eg.pS]) 
to the order in q indicated. 



B. Definition of the height profile 

An important subtlety in the preceding anal- 
ysis is the fact that the location of the interface, 
i.e. the value of the height function ft.(r||), can- 
not be defined unambiguousl]^. A certain pro- 
cedure must always be formulated to determine 
ft.(r||). It turns out that the choice for h{r^^) in- 
fluences the density profile pi{z) which, in turn, 
determines the value of the bending parameters 
k and £k- 

We explicitly consider two canonical choices 
for the determination of /i(^| ); the crossing con- 
straint (cc) and the integral constraint (ic) ^^i^^ . 
Other choices are certainly possible and equally 
legitimate as long as they lead to a location of 
the dividing surface that is 'sensibly coincident' 
with the interfacial region?!. In this context we 
like to mention the work by Tarazona et al^, 
who propose a 'state of the art' manner to de- 
fine the location of the interface based on the 
distribution of molecules rather than the molec- 
ular density alone. 
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In the crossing constraint, h{f\^) is defined as 
the height where the fluctuating density equals 
some fixed value of the density that lies in be- 
tween the limiting bulk densities, say p(r) = 
po{z = 0): 



p(f||,z = /i(r||)) = po(0). (cc) 



(16) 



Using this condition in Eq. ^ , one finds the fol- 
lowing constraint for pi(z) — p1'^{z) 



pr(o) = o. 



(17) 



In the integral constraint, h{r\\) is defined by the 
integral over the fluctuating density^^ 



hifn 



1 



dz[p{r) - Pstcpiz)] ■ (ic) (18) 



The quantity studied in experiments and sim- 
ulations is the (surface) density-density correla- 
tion function. It is an integral into the bulk re- 
gion to a certain depth L of the density-density 
correlation function: 



1 



L L 

dzi / dz2 



(23) 



(Ap)2 

< [P(^l) - Pstcvizi)] [P{r2) - Pstop(22)] > ■ 

When we insert the general expression for p{'r) 
as given by Eq.Q into Eq. ([23| . one finds that 

S{r\\) - <M^1,||)M^\||)> (24) 

oo 

-■^ Jdzpi{z) <h{fi^ll)Ah{r2,\\)> 



With this condition inserted into Eq.®, one 
now finds that pi{z) — pf{z) is subject to the 
following constraint 



dzpi'^iz) = 0. 



(19) 



We show in Section 4 that the ambiguity in lo- 
cating the dividing surface translates into the 
density profile pi{z) being determined up to an 
additive factor proportional to Pq{z)'^'^. In par- 
ticular, pf{z) and Pi^{z) are related by 



pr{z) = pT{z) + ap',{z)^ 



(20) 



The value of the constant a can be determined 
by integrating both sides of the above equation 
over z 



oo 



(21) 



where we can neglect a term < Ah Ah > to 
the order in the curvature expansion considered. 
Furthermore, we have assumed that L is suffi- 
ciently large so that we can approximate 



dz p'q{z) 



-L 
L 



dz pi{z) 



-Ap, 



dzpi{z) 



(25) 



Rather than 5(r||), we consider its Fourier 
Transform, S{q), which we shall term the sur- 
face structure factor. 



(26) 



Siq) - /dr||e-''-'-|i5(r||) 



S'/t/i(g) + I dz pi{z) Shh{q) 



One may further show that the ambiguity in 
the determination of pi{z) is of influence to the 
value of the bending parameters k and ik- In 
Section 4 we show that because p^i{z) and p'i'^{z) 
are related by Ea. ipD]) . we have for the bending 



parameters'^ 



1 ic 1 cc _ 

k = k — acr , 

fc. in(£r) = ks Hei^) - 



(22) 



Naturally, all experimentally measurable quan- 
tities cannot depend on the choice made for the 
location of the height function h{f\\). The im- 
plication is that it is necessary to formulate pre- 
cisely the quantity that is determined experi- 
mentally and verify that its value is independent 
of the choice for h{f\\). This is explicitly shown 
next. 



We now verify that S{q) is independent of the 
choice for /i(r||) by determining S{q) using both 
the integral constraint and crossing constraint. 
For simplicity, we consider the case of short- 
ranged forces only (the verification for the case 
of long-ranged forces follows analogously) . The 
surface structure factor using both constraints 
is given by 



5-(g) = 



k^T 



a k^T q^ 



a q^ - 



V k^" q^ + 
k^Tk" 
&^ 

knT 



erg- 4 
a k^T 



q- 



a q^ ~ 
kBT 



k^T¥ 



0{q^), 



Oiq'), 



(27) 
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FIG. 1: Sketch of the fluctuating density proflle as 
a function of z; the height h = h[r\\) is the distance 
over which the intrinsic density proflle po{z) (dashed 
hne) is shifted. 



where we have used the explicit expression for 
Shh{q) in Eq.((T2l) together with Eq.dlSl). On 
account of the fact that fc*^ = k'^'^ — aa, one finds 
that S'"'{q) = S''''[(i) = S{q) as required. 

This analysis shows that S'(q) equals the 
height-height correlation function when the in- 
tegral constraint is used to define the location of 
the height profile, i.e. 

S{q)^S^M. (28) 

It is therefore convenient, but by no means nec- 
essary, to use the integral constraint to define 
the location of the dividing surface. 

Finally, we consider the contribution of "bulk- 
like" fluctuations to the fluctuating density pro- 
file which are predominantly present at short 
wavelengths, qd>\. 



C. Bulk-like fluctuations 

Adding short wavelength, bulk-like fluctua- 
tions to the fluctuating density, the full picture 
that emerges for p(f) is that schematically de- 
picted in Figure [2 It can be described as: 

p{v) = p^{z)~p'^{z)h{r{)-P^ ^h{r\{)^ip,{r) 

(29) 

where Spb{r) represents the bulk- like fluctua- 
tions. We shall consider only small fluctuations 
so that < 6pb >= and assume that there are 
no correlations between height fluctuations and 
bulk- like fluctuations, < hSpt >= 0. When we 
insert the expression for p(r) as given by Eq. (P5)) 
into the expression for S{r^\) in Eq. ((23)) . one 
finds that 

5(r||)= Slli^) (30) 



L oc 

~''(Ap)2 /^^^ jdzi2 <5pb{ri)5ph{r2)> ■ 

— L — oo 

Here we have made a further approximation by 
replacing the integration over 22 from —L to L 
by an integral over 212 from —00 to 00. The in- 
tegral over zi that is left gives rise to a term that 
increases linearly with L. That means that the 
bulk- like contributions to S'(rii) eventually dom- 
inate the height fluctuations when L becomes 
larger. To study surface fluctuations via S{r\\) 
it is therefore important that on the one hand 
L is sufficiently large in order to make the ap- 
proximations in, e.g., Eq. (|25p but on the other 
hand not so large as to completely dominate the 
contribution from surface height fluctuations. In 
the next section we show how these two condi- 
tions pan out for the circumstances under which 
the simulation results are obtained. 

A further issue is that the bulk density corre- 
lation function < Spb Spb > differs in either phase 
(liquid or vapor). When one then considers the 
integral over zi , it seems appropriate to approxi- 
mate < Spb Spb> by the density correlation func- 
tion in the bulk liquid region: 

<Spb{ri)Spb{f2)>^ pj [ge{r) - l]+piS{fi2) , 

(31) 

and introduce an L-dependent prefactor JVl to 
account for the integral over zi. The surface 
structure factor thus becomes 

Siq)^Slliq)+N-LSb{q), (32) 

with the bulk structure factor Sb{q) defined as 

Sb{q)^l + p, ydri2e-*«"- . (33) 

This approximation may be justified by arguing 
that close to the critical point there is no distinc- 
tion between the two bulk correlation functions, 
whereas far from the critical point the contribu- 
tion from the bulk vapor can be neglected since 
Pv~0- 

The value for the L-dependent prefactor A/l 
may be determined from a fit to the limiting be- 
havior of S{q) at q — >oo. For an explicit evalua- 
tion of Sb{q), we have taken for gi{r) the Percus- 
Yevick solution^^ for the hard-sphere correlation 
function, gi,{r)=glJ [r] pi). 

III. COMPARISON WITH MONTE 
CARLO SIMULATIONS 

In this section, the surface structure factor in 
Eq. ([5^ is compared to results from Monte Carlo 
simulations by Vink et al.^. The system con- 
sidered consists of a mixture of colloidal parti- 
cles with diameter d and polymer particles with 
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FIG. 2: MC results by Vink (Ref. i) for the surface 
structure factor (in units of d*) versus q (in units 
of l/d) for various values of the integration limit 
L/W = 1, 2, 3, 4. The dashed line is the capillary 
wave model. In this example e= 1.8, rjp — 1.0, and 
the colloidal particles are used to define zq- 



diameter 2i?g. The colloid-colloid and colloid- 
polymer interactions are considered to be hard- 
sphere like, whereas polymer-polymer interac- 
tions are taken ideal. The presence of poly- 
mer induces a depletion attraction between the 
colloidal particles which may ultimately lead to 
phase separatiorti^ii^iilii^. The resulting inter- 
face of the demixed colloid-polymer system is 
studied by Vink et al.^ for a number of polymer 
concentrations and for a polymer-colloid size ra- 
tio parameter e=l -I- 2Rg/d = 1.8. 

To study the interfacial fluctuations, Vink et 
al. introduce the local interface position as^: 



1^ 

-L 



(34) 



where p(f) can be taken to be either the colloid 
or polymer density. The integration limits ±L 
are inside the bulk regions, but different values 
for it are systematically considered"^. One may 
easily verify that the correlations of the local in- 
terface position are exactly described by the sur- 
face structure factor defined earlier in Eq. 



<ZG{ri^\\) ZGif2.\\)>= 5(r||). 



(35) 



In Figure[21 typical results for the Fourier trans- 
form of the surface structure obtained in the 
MC simulations of Vink are shown (Figure 13 
of ref. 3). In this example the integration limit 
is varied, L/W = 1, 2, 3, 4, where W is some 
measure of the interfacial thickness. One clearly 
observes that when L/W is too small, the re- 
sults do not match the classical capillary wave 
behavior for small q (dashed line), and that the 
contribution from bulk-like fluctuations at high 
q increases with L/W. 



FIG. 3: MC results by Vink et al. (Ref. ;3;) for the 
surface structure factor (in units of d'^) versus q (in 
units of 1/d). The dotted line is the capillary wave 
model, the dashed line is the combination of the 
capillary wave model and the bulk correlation func- 
tion, and the drawn line is the combination of the 
extended capillary wave model and the bulk corre- 
lation function. In this example e = 1.8, rjp = 1.0, 
L/W = 3, and the colloidal particles are used to de- 
fine zq. 



In Figure [31 we consider the result from Fig- 
ure [2] for L/W = 3. For small q the results 
asymptotically approach the result of the classi- 
cal capillary wave model (dotted line) with the 
value of a taken from separate simulations. The 
dashed line is the combination of the capillary 
wave model with the bulk correlation function: 



S{q)^ 



aq^ 



-ULSbiq)^ 



(36) 



The value of Ml is chosen such that it matches 
the g— s-oo limit for S{q) in Figure [3] One finds 
that Eq. ([55)) already matches the simulation re- 
sults quite accurately except at intermediate val- 
ues of q, qd~ 1. 

As a next step, we investigate whether the in- 
clusion of a bending rigidity is able to describe 
the simulation results at these intermediate val- 



ues: 



S{q) 



kBT 



a q-^ 



+ MLSb{q). (37) 



The bending rigidity describes the leading order 
correction to the classical capillary wave model 
in an expansion in q^. Its value is therefore ob- 
tained from analyzing the behavior of S{q) when 
qd < 1. The fact that the simulation results 
in Figure [3] are systematically above the cap- 
illary wave prediction in this region, indicates 
that the bending rigidity thus obtained is nega- 
tive, fc**^ < 0. Unfortunately, a negative bending 
rigidity prohibits the use of Eg. ([57)) to fit the 
simulation results in the entire g-range since the 
denominator becomes zero at a certain value of 
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q. It is therefore convenient to rewrite the ex- 
pansion m 
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in Ea. ((37|) in the foUowing form: 



S(q)^^il-^q^ + ...)+ULS,iq), (38) S(q) 
aq^ a . 



which is equivalent to Eg. ([57)1 to the order in q 
considered, but which has the advantage of be- 
ing well-behaved in the entire g-range. Other 
forms to regulate S{q), that are equivalent to 
Eq. ([37l) to the order in q considered, may cer- 
tainly be formulated. In analogy with a simi- 
lar treatment of capillary waves by Parry and 
coworkers^i in the context of wetting transi- 
tions, one might suggest that the appearance of 
a negative bending rigidity indicates the miss- 
ing of a correlation length that would replace 
Eq. ((38)) with an explicit formula valid for all val- 
ues of q, not just to the order in q considered. 

The above form for S{q) in Eq. pS)) . with the 
bending rigidity used as an adjustable parame- 
ter {k= - 0.045 k^T), is plotted in Figure [3] as 
the drawn line. Exceptionally good agreement 
with the Monte Carlo simulations is obtained. 
In Table 1, we list values of the bending rigidity 
obtained for a number of polymer volume frac- 
tions, rjp. These values are the results of fits 
of S{q) from Monte Carlo simulations for sev- 
eral system sizes and for several values oi L/W, 
with the error estimated from the standard de- 
viation of the various results. For the L/W = 1 
and L/W = 2 curves (see Figure [2]), one needs 
to adjust for the fact that the capillary wave 
limit is not correctly approached at low q. For 
the results in Figure [H one ultimately obtains 
for the bending rigidity k = - 0.040, - 0.040, - 
0.045, - 0.060 fceT, for L/W^l, 2, 3, and 4, 
respectively. 

In Table 1, it should be reminded that, rather 
than the true polymer volume fraction in either 
phase, rjp should be interpreted as the polymer 
volume fraction of a reservoir fixing the polymer 
chemical potential^. Furthermore, the "liquid" 
is defined as the phase relatively rich in colloids 
and the "vapor" as the phase relatively poor in 
colloids. 

The excellent agreement between Eq. ([55)1 and 
the MC simulations is even more apparent in 
Figure [4] where the results in Figure [3] are re- 
drawn on a linear scale. In Figure[3]we also show 
the simulation results^ and the corresponding fit 
using the polymer particles to define the location 
of the interface. As the polymer-polymer inter- 
actions are considered ideal, the bulk structure 
factor Sb{q) = l in this case. 

It is important to note that, effectively, the 
inclusion of a bending rigidity in the capillary 
wave model results in the presence of an addi- 
tive factor in S{q) that is adjusted, see Eq. ([38|) . 
The determination of the value for fc'^ from the 
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FIG. 4: MC results by Vink et al. (Ref. i) for 
the surface structure factor (in units of d**) versus 
q (in units of 1/d) using the colloidal particles (cir- 
cles) and polymer particles (triangles) to define zq- 
The dotted line is the capillary wave model and the 
drawn lines are the combination of the extended cap- 
illary wave model and the bulk correlation function. 
In this example e= 1.8, rjp = 1.0 and L/W= 3. 



Vp 


Vi 


r)v 




k 




0.9 


0.2970 


0.0141 


0.1532 


-0.045 (15) 


0.40 


1.0 


0.3271 


0.0062 


0.2848 


-0.07 (2) 


0.50 


1.1 


0.3485 


0.0030 


0.4194 


-0.10 (3) 


0.49 


1.2 


0.3647 


0.0018 


0.5555 


-0.14 (3) 


0.50 



TABLE L Listed are the simulation results (Ref. 0) 
for the polymer volume fraction rjp, liquid and vapor 
colloidal volume fractions, rji and surface tension 
a (in units of ksT/d?), bending rigidity k (in units 
of fceT; in parenthesis the estimated error in the last 
digit), and — k/a (in units of d). 



behavior of S{q) near q — Q, therefore requires 
one to take into account the presence of the 
bulk-like fluctuations since they also contribute 
as an additive constant, J\fLSb{0), near q = 0. 
This means that even though the MC simulation 
results of Vink et al^ are very accurately de- 
scribed by Eq. , the resulting value obtained 
for fc*"^ sensitively depends on the theoretical ex- 
pression used for <S'f,(0). Here we have simply 
approximated the bulk correlation function by 
the Percus-Yevick hard-sphere expression in the 
liquid^^, but one could imagine more sophisti- 
cated expressions leading to a somewhat differ- 
ent value for k^'^. 

In the next sections we investigate whether 
the values for the bending rigidity obtained from 
the simulations (Table 1), can also be described 
in the context of a molecular theory. 
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IV. DENSITY FUNCTIONAL THEORY 

Our task in this section is straightforward. 
Using the expression for p[r) given in Eq.Q, 
we determine AO and the resulting a{q). To 
achieve this, we need a model for the free energy 
and define a procedure to determine the density 
profiles, Pq{z) and pi{z), that are present in the 
expression for p{f). We choose to perform these 
tasks in the context of density functional theory 
(DFT). 

In the density functional theory for an inho- 
mogeneous system that we consider— i^iiii^, the 
free energy is given by the free energy of the 
reference hard sphere system augmented by an 
integral, non-local term that considers the at- 
tractive part of the interaction potential, U{r) = 
[/hs(r) + t/att(r). 



(39) 



dfi dr2 U^tt{r) p{ri)p{r2) 



For explicit calculations, ghsip) taken to be of 
the Carnahan-Starling form^^: 

5hs(p) = ksTp \n{p) + ^bT'p^^y— - pp, 

(40) 

where rj = [tt / 6) p . In the uniform bulk region, 
the free energy equals 



V 



gip) = ghsip) - ap^ , 



(41) 



with the van der Waals parameter a given byi 
1 



df 12 Ui,ttir) . 



(42) 



The integration over ri2 is restricted to the 
region r > d. This is not explicitly indi- 
cated; instead, we adhere to the convention that 
the attractive part of the interaction potential 
Uutt ir) — when r <d. The chemical potential 
p is fixed by the condition of two-phase coex- 
istence, p = Pcocx, which implies that Pcocx, Pvi 
and Pi are determined from the set of equations: 
g'iPv)^0, g'iPe)^0, and giPv) = giPe)^-P- 

To determine the change in free energy due 
to density fluctuations, we insert the expression 
for pir) given by Eq.® into the expression for 
n in Eq.dSS]). One then finds for A0 = - crA 



AO = i Jdn {gLiPo)Piizi)' [Ahiri,\\)y 

+ ^ Jdri jdf 12 U^ttir) x 
|-2po(zi)po(22) [^(7^2,11) - hiri^\\)Y 



+4pi(^i)Po(^2) A/i(fi,||) [hir2,\\) ~ Kru])] 
+Piizi)piiz2) A/i(ri^||)A/i(r2^||)} . (43) 

Even though the derivation is somewhat differ- 
ent, this expression equals that given by Mecke 
and DietricblS apart from a gravity term that 
was included in their expression. To cast AO in 
the form of Eg.pT]). we take the Fourier Trans- 
form. One then finds for cr(g)^ 



<l) 



dzi / dzi2 



'iq, Z12) -^0(^12) 



— 00 —00 



X [Poizi)Poiz2) ~ q'^ Piizi)p'oiz2) _ 

00 00 

„2 



+ ^ I dzi I dzi2 uiiq^zu) piizi)piiz2) 



+ X /^^i 5hs(Po)pi(2i)^ 



(44) 



Here we have defined the (parallel) Fourier 
Transform of the interaction potential 



u;iq,zi2) = /drii e-^'-'-ii C/att(r) (45) 



= 2iT jdr\\r\\Jniqr\\)Ua^ttir) . 


As a first step, we determine the leading con- 
tribution to aiq) given by the surface tension of 
the planar interface, a = aiq = 0). Then, one 
needs to consider the two leading contributions 
in the expansion of cj(q, Z12) in q^: 

w(g, Z12) = uJoizi2) + uj2izi2)q^ + . . . , (46) 
where 



^0(^^12) = fdf\\ Ua,nir) 



0.2(^12) ^ -J dr\\ r^C/att(r). (47) 



The surface tension thus becomes 



dzi dzi2 uj2izi2) p'oizi)p'Qiz2) ■ (48) 



-00 —00 



The (planar) density profile poiz), featured 
in the above expression for ct, is determined 
from minimizing the free energy functional Q[p\ 
in Eq. (j39p in planar symmetry. The Euler- 
Lagrange equation that minimizes 0[p] is then 
given by: 



g'hsiPo) - J dfi2 U^ir) Poiz2) 

00 

= - dzi2 ujoizi2) poiz2) , (49) 
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which can be solved exphcitly to obtain pq{z) 
and thus a. 

The evaluation of further contributions to 
(T{q) requires one to determine the density pro- 
file pi{z). Just like Pq{z), one would like to de- 
termine the density profile pi{z) from a mini- 
mization procedure. One then has to determine 
the energetically most favorable density profile 
for a given curvature of the surface^^. This turns 
out to be not so straightforward, since one then 
has to specify in what way the curvature is set 
to its given value. Several approaches have been 
suggested, which we shall now discuss. 

• Mecke and Dietrich approach. In this 
approach a certain form for pi{z) is di- 
rectly hypothesizedi^: 



ApUniz/O^ 



(50) 



with ^ the bulk correlation length and 
/H(a;) = a;sinh(a;/2)/cosh^(x/2). The co- 
efficient Ch in this expression can be used 
as a fit parameter. This practical ap- 
proach is certainly legitimate, but one 
would like to also be able to formulate a 
molecular basis for this expression. 

• EquiUbrium approach. Rather than 
the surface being curved by surface fluc- 
tuations, in this approach the interface is 
curved by changing the value of the chem- 
ical potential to a value off- coexistence. 
One then considers the density profile of 
a spherically or cylindrically shaped liquid 
droplet in metastable equilibrium with a 
bulk vapor^"'^^. This approach is equiva- 
lent to adding an external field to the free 
energy 

n'[p] = n[p]+ [dfV,^t{r^p{r) 



(51) 



= ri[p] - JdrApp{f) 



where lS.p = p — pcocy.- The downside of the 
equilibrium approach is that the external 
field Vcxt(?') = — A/u is uniform through- 
out the system and thus also affects the 
bulk densities far from the interfacial re- 
gion. This seems inappropriate for the de- 
scription of the density fluctuations con- 
sidered here since we have that p — Pcocx 
and the bulk densities are unaltered by the 
curvature of the surface fluctuations. 

• Local external field. In this approach, 
one again adds to the free energy an ex- 
ternal field, but, to ensure that the bulk 
regions are unaffected, one assumes that 
it is peaked infinitely sharply at z = O^li^: 



In this case, the external field only acts 
as a Lagrange multiplier in the minimiza- 
tion procedure to ensure that the curva- 
ture A/i(r||) is set to a certain value; it is 
not included in the expression for the free 
energy. The downside of this method is 
that the resulting density profile pi{z) has 
a discontinuous first derivative at z = 0, 
which is, from a physical point of view, not 
so appealing--. Furthermore, the discon- 
tinuous nature of pi{z) prohibits an ana- 
lytical simplification using a gradient ex- 
pansion. 

In the present approach, we suggest to add an 
external field acting as a Lagrange multiplier 
that is unequal to zero only in the interfacial 
region (the bulk densities are unaffected), but 
which is not infinitely sharp-peaked. It seems 
natural to choose a peak-width of the order of 
the thickness of the interfacial region. It thus 
seems convenient to choose Vcnti^ oa p'q{z): 



Kxt(r)-Ap;,(z)A/i(f||). 



(53) 



This choice for Vcxt(^ constitutes our funda- 
mental 'Ansatz' for the determination of pi{z). 
The Lagrange multiplier A is not a free parame- 
ter but set by the imposed curvature, as demon- 
strated below. 

The addition of an external field to the free 
energy results in the following Euler-Lagrange 
equation: 

g'M = -Jdr 12 UMr) P{r2) - Kxt(r) . (54) 



Using the external field given in Eq. (j53p . we in- 
sert the fluctuating density given by Eq. ^ into 
the above Euler-Lagrange equation. In order 
for the resulting equation to hold independently 
of the value of h{r\\) or A/i(r||), one finds, be- 
sides Eq. (P^ . the following equation to deter- 
mine pi{z): 



5hs(Po)/Oi(2i) = - ldzi2 ^0(^12) Pi (22) 



+2 dzi2 ^^2(2:12)^0(22) + 2 Ap^(zi) . (55) 



The value of the Lagrange multiplier can be de- 
termined by multiplying both sides of the above 
expression by /5o(zi) and integrating over zi. 



A 



dz [p'^{z)f 



(56) 



K=xt(rl = A5(z)A/i(n|) 



(52) 



One may now verify that if pi{z) is a particular 
solution of Eq. ((55|) that then also pi{z) -\- a p'q{z) 
is a solution on account of Eq. . 
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It is convenient to use the Euler-Lagrange 
equation in Eg. ([55]) to remove the explicit ap- 
pearance of g{[s{po) in the expression for a{q) 
in Eg. p^ . The resulting a{q) is written as the 
sum of a term that depends only on the density 
profile po{z) and one term that also depends on 
the density profile pi{z) 



a{q) = aoiq) + h q^ + 0{q^) 



(57) 



with 



h = - 



dzi / dzi2 



U}{q, Z12) - ^0(2:12) 



xPo (^1)^0(2^2), 

dzi I dzi2 UJ2{zi2) Pi{zi)p'q{z2) 

(58) 



+ 2 I dzi pi{zi)p'o{zi) . 



With the above expression for ki it is now also 
possible to verify that when the density profile 
pi{z) is shifted by a factor Q!Pq(z), that the re- 
sulting effect on the bending parameters is that 
given by Eg. ([22]) . 

The procedure to determine <7{q), and there- 
fore (7, fc, fcs, and £k, is now as follows: as- 
suming a certain form for the attractive part 
of the interaction potential, po{z) is obtained 
from solving Eq. ipS)) . which is then inserted into 
Eg. ((55)) to solve for pi{z) exphcitly. The two 
density profiles thus obtained are inserted into 
Eg. ([58]) to yield cro{q) and fci. This procedure 
is carried out in the next two sections consid- 
ering short-ranged forces and long-ranged forces 
(C/(r) cx l/r^)- In general the density profiles 
Po{z) and pi{z) need to be determined numer- 
ically. We shall, however, also provide an ap- 
proximation scheme, based on the gradient ex- 
pansion, that is exact near the critical point, but 
which also gives an excellent approximation far 
from it. 



where m is the van der Waals squared-gradient 
coefficient^ 

m = Jdri2r^UMr) (61) 

00 00 

= ^ ^2(212) = ^ Jdzi2 zf2UjQ{zi2) . 



In the gradient expansion, the Euler-Lagrange 
equation in Eg. (|55p for pi{z) reduces to 



m—p——pi{zi) = mpi(zi) 
Po(^i) 



(62) 



+ /rf2:i2 ^2(2:12) Po(z2) + A Po(zi), 



where we have used Eg. ([SO)) to replace 5(,'g(po)- 
First, we consider the determination of the 
density profile po{z). The gradient expansion 
becomes exact near the critical point where g{p) 
takes on the usual double-well form 



l{p)+P = 



{Ap)H' 



[p- p,f{p- p,f . (63) 



Using this form for g(/o), the solution of the 
Euler-Lagrange equation in Eg. ([50]) gives the 
usual tanh-form for ^0(2)^: 



Po{z) = + Pv) 



Ap 



tanh(z/20 , (64) 



with the bulk correlation length ^ a measure of 
the interfacial thickness. 

Even though the tanh-form for the density 
profile po{z) is derived assuming proximity to 
the critical point, it turns out that it also pro- 
vides a good approximation away from it when 
one determines the value of ^ by fitting the sur- 
face tension to its form near the critical point. 
In the squared-gradient approximation, Eq. 
reduces to the familiar expression^ 



a = 2 m dz Pf){z) 



A. Gradient expansion 

The gradient approximation^ is based on the 
assumption that the spatial variation of the den- 
sity profile is small, i.e. 



piz2) = p(zi)+zi2 p'izi) + ^ p"{zi) + . . . (59) 

In the gradient expansion, the Euler-Lagrange 
equation in Eg. (|49p for pq{z) reduces to 



g'{po) = 2mp'd{z), 



(60) 



= J dp V9{P)+P- (65) 

On the one hand, the surface tension can be 
determined from the above approximation using 
the full form for g{p) given in Eqs.([40l) and (|4ip : 

Pe 

(7 = 2 \/m jdp \k^Tp\n{p) 



(66) 



-l-fceT p _ Mcoox/O - ap^ + p] ^ 
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On the other hand, near the critical point g{p) 
takes on the double- weU form in Eq. ((63|) and a 
is calculated as 



(67) 



Now, we define the value of ^ such that the two 
expressions for the surface tension in Eqs.(|66p 
and (|67p are equal. This gives for ^: 



■m{ApY 

3cr 



(68) 



with a given by Ea. (|66|) . 

Next, we turn to the evaluation of pi{z) from 
Eq. ([62|) . This requires one to make a dis- 
tinction between short-ranged forces and long- 
ranged forces. 
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FIG. 5: Surface tension in units of ksT/d^ versus 
the volume fraction difference, Arj. In this example 
e = 1.8; symbols are numerical results, the drawn 
line is the gradient expansion approximation, filled 
symbols are results from the MC simulations by 
Vink et al. (Ref.i). 



V. DFT: SHORT-RANGED 
INTERACTIONS 

Although the analysis below is quite generally 
valid for all short-ranged interaction potentials, 
whenever we show explicit results, we consider 
for Uattir) the Asakura-Oosawa-Vrij depletion 
interaction potential Udcpi'r) as an example^: 



2 (£-1)3 



(69) 

where the intermolecular distance is in the range 
l<r/d<e. Interaction parameters based on the 
depletion potential are listed in the Appendix. 
The strength of the depletion interaction po- 
tential as determined by the polymer volume 
fraction rjp determines the location in the phase 
diagra m^^'^? ; for comparison with other results, 
it is, however, more convenient to use the col- 
loidal density difference A77 = 77^ — 77^ as thermo- 
dynamic variable^. 

In Figure O the surface tension is shown as a 
function of A77. The open circles are obtained 
from numerically solving the Euler-Lagrange 
equation in Eg. ([49)1 for po{z) and inserting the 
result into Eq . ((48|) . The drawn line is the gra- 
dient expansion approximation for a in Eq. (|66p . 
Also shown are results from the MC simulations 
by Vink et al^. The gradient expansion gives 
a very good approximation to the numerical re- 
sults and is in good agreement with the simula- 
tions. 

For short-ranged forces the expansion in of 
the expression for uj{q, Z12) as defined by Eq. (|45ll 
can be continued to 0{q'^): 



ui{q,zi2) = Wo(2i2)+W2(zi2) 9^+^^4(212) 9"* + - 



where 



1^4(212) 



1 

64 



dr\\rlU^tt{r). (71) 



With this expansion, a{q) in Eq. ([58|) can now 
be written in the form of Eq. (fT3|) 

a{q) = (T + kq^ + 0{q^) 

= <J + koq^ + hq^ + 0{q^), (72) 

with the bending rigidity k — k^^-ki and 



00 00 



fco 



fcl 



dzi dzi2 UJi{zi2) Po(2l)/0o(2;2) : 



00 00 



dzi dzi2 UJ2{Z12) Pi{zi)Po{z2) 



-00 —00 



dzi pi{zi)pq{zi) . 



(73) 



Next, we proceed to evaluate these expressions 
in the gradient expansion. 



A. Gradient expansion for short-ranged 
forces 



In the gradient expansion, ka and fci in 
Eq.dlSl) reduce to: 



fco = -— jdz p'Q{zf 



-00 



ki = —2m 



(70) 



Jdzpi{z)p'o{z), (74) 



12 



where we have used the fact that to leading order 
in the gradient expansion A « —2m, and where 
we have defined'^^ 



B 



1 

60 



df 12 r Ui,tt{r) 



(75) 



oc oo 
= -2 Jdzi2 UJ4{zi2) = ^ Jdzi2 zl^U02[zi2). 

— oo — oo 

Inserting the tanh-form for po{z) into Eg. ((74|) . 
one directly obtains for fco 



fco = - 



12 e 



(76) 

4 m 



where we have used the expression for ^ in 
Ea. ((68| to rewrite as the latter expression. 

To evaluate ki in Eg. ([7^ . we need to deter- 
mine pi{z) from the Euler-Lagrange equation 
in Eq. ([S^ . For short-ranged forces, Eg. lp^ re- 
duces to 



TO— ^Pl(z) 



mp'i{z) + Bpl'{z) + pBp'^{z), 

(77) 



where we have defined 

Sdz p'^{zf 



/3: 



jdz p',{zf 



(78) 



Using the tanh-profile for po{z) in Eg . ([64| . one 
has /9=l/(5^^) and finds for pi{z) from solving 
the differential eguation in Eg. ((77)) : 



3B Ap ln(cosh(z/20) 



10 m ^ cosh^(z/20 



(79) 



The above profile corresponds to that obtained 
using the crossing constraint. The profile corre- 
sponding to the integral constraint follows from 
Pi' = PT + Po (Eq-dSOl)), with a determined by 
Eg. (EI]). This gives 



pTi^) 



3B Ap [1 - ln(2 cosh(z/2^))] 



10 m ^ 



(80) 



cosh^(z/20 

In Figure [51 typical volume fraction profiles 
rji (z) — (tt/G) d^ pi (z) are shown for the crossing 
constraint and the integral constraint. The sym- 
bols are the profiles obtained from numerically 
solving the Euler-Lagrange eguation in Ea.(|551). 
whereas the drawn lines are the approximate 
profiles in Eqs. ([79)) and ([80)1 obtained from the 
gradient expansion. 

Inserting the density profiles p^i{z) and pf^iz) 
into the expression for ki in Eg. ([7^ . one obtains 



ki" 



B {Apf 



ln(2) 



B<j 



BjApf 



1 - - ln(2) 

o 

Ba 
5 m 



(81) 
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FIG. 6: Curvature correction to the volume frac- 
tion profile 771(2) as a function of z (in units of d) 
using the integral constraint (circles) and crossing 
constraint (squares). In this example At] = 0.25, 
e = 1.8; symbols are numerical results, the drawn 
lines are the analytical profiles from the gradient 
expansion. 




FIG. 7: Contributions to the bending rigidity fco 
and fci in units of k^T versus the volume fraction 
difference, A77. In this example e = 1.8; symbols are 
numerical results, the drawn lines are the gradient 
expansion approximation. 



In Figure [3 feg, k^f and kf^ are shown as a 
function of Arj. The open symbols are obtained 
from numerically solving Egs. iPS]) and ([55]) to 
obtain pq{z) and pi{z) and inserting the result 
into Eg. ([73)) . The drawn lines are the gradi- 
ent expansion approximation for ko in Eg. ()76[) 
and ki in Eg. ([8T)) . Adding the results for ki 
in Eg. ([8T)) to fco in Eg. ([76)1 . one obtains for the 
bending rigidities 



B{Ap) 



2 r 



5 

12 



ln(2) 



Ba 



^ - ^ ln(2) 



B{Apf 
60^ 



Ba 
'20 m 



(82) 
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FIG. 8: Bending rigidity in units of ksT versus 
the volume fraction difference, Arj, using the in- 
tegral constraint (circles) and crossing constraint 
(squares). In this example e = 1.8; open symbols 
are numerical results, the drawn lines are the gra- 
dient expansion approximation, filled circles are the 
MC resuhs by Vink et al. (Ref.i); the dashed line 
is the fit yj-k/a ^ 0.47 d. 



In Figure [8l the bending rigidity is shown as 
a function of A77. The open symbols are the 
numerical results. The drawn lines are the 
gradient expansion approximations in Eq. (|82p . 
y^FV^ « 0.378 d and sj-k^'^jo « 0.131 d. 

The resulting values for both bending rigidi- 
ties, k'^'^ and fc*'^, are negative in line with the 
simulation results of Vink (filled circles). How- 
ever, the value of fc'^, which is the relevant value 
when we compare with the simulations, is sig- 
nificantly less negative. As stressed earlier, the 
bending rigidity depends on the constraint used 
to define the height profile h{r\\) through the 
contribution to k coming from k\. We demon- 
strated that in the crossing constraint k'^ is neg- 
ative whereas in the integral constraint k"^-^ is 
positive. One could very well imagine that a dif- 
ferent constraint used to determine the height 
profile h(r\\ ) might lead to a bending rigidity 
that is positive'*. For the integral constraint the 
two contributions to k from fco and k'^{' nearly 
cancel leading to a value for fc*'^ which is barely 
negative. Unfortunately, this makes the value of 
fc*"^ sensitively dependent on the precise model 
used to determine p\{z)- 

An important point concerns the scaling be- 
havior of the bending rigidity. The expressions 
in Eg. ([5^ indicate that the bending rigidity 
vanishes near the critical point with the same 
exponent as the surface tension, i.e. 



k oc cx era , 

m 



(83) 



Note that for the depletion potential, the ratio 
B/m only depends on the size ratio parameter 
e {B/m « 0.342 for e= 1.8) but is indepen- 



dent of rjp (or A77), see the Appendix. Both 
contributions to the bending rigidity, fco and ki , 
show the above scaling behavior and both should 
therefore be taken into account. 

The scaling result in Eq. ([55|l should be con- 
trasted with the usual assumption that k(xa^^, 
i.e. the bending rigidity approaches a finite, 
non-zero limit at the critical point. This scaling 
behavior is, for instance, obtained for the bend- 
ing rigidity fceq determined from analysing the 
surface tension of a spherically or cylindrically 
shaped liquid droplet in metastable equilibrium 
with a bulk vapor— i^. In the gradient expan- 
sion, one has^S: 

fceq = -i(vr2-3)m(Ap)2e. (84) 

It is perhaps important to discuss more broadly 
this result in the context of previous work on 
the virial approach^! to the bending rigidity 
(and other curvature parameters). The virial 
expression for the bending rigidity is generally 
valid, but it is important to realise that it fea- 
tures the way in which the pair density depends 
on curvature^'*. When a mean-field, squared- 
gradient approximation is subsequently made^, 
this translates into the expression for the bend- 
ing rigidity to depend on the way in which the 
density depends on curvature, i.e. it features the 
profile pi{z). Therefore, even though the expres- 
sions for the bending rigidity in the equilibrium 
approach and the fluctuating interface approach 
are the same, they might lead to different values 
(and scaling behavior) of the bending rigidity 
due to the fact that the density profiles pi{z) 
are different in these two cases. 

It is interesting to also compare with the ap- 
proach by Mecke and Dietrich^^. Even though 
the goal in ref. [l9; is to consider long-ranged 
forces, one may use the expression in Ea. (|50p for 
pY^{z) inserted into Eg. ([43)1 to determine the 
leading correction to cr(q) also for short-ranged 
forces. The gradient expansion then gives: 

^" (7r2-|-16)lm(Ap)2^. 



kuB 



6 24 7r2 J 

(85) 

The prefactor is negative (as long as < Ch < 
1.526) in line with the results obtained here. 
Again, the scaling behavior - equal to that of 
fccq ~ is essentially different than our prediction 



in Eq. 

Finally, we like to mention an expression for 
the bending rigidity that is derived from the gen- 
erally valid virial expression^!, in which the as- 
sumption is made that the width of the inter- 
facial profile is much smaller than the molec- 
ular diameter d ^ (^24^33 ^ r^^ie implication is 
that po{z) = pstcp{z), the sharp-profile approx- 
imation, and pi{z) = 0. The sharp-profile ex- 
pressions for the surface tension (also known 
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FIG. 9; Bending rigidity in units of ksT versus the 
volume fraction difference, Ary. In this example e — 
1.8; filled circles are the MC results by Vink et al. 
(Ref.i), the drawn line is the virial expression with 
the sharp-profile approximation, the dashed line is 
the equilibrium result, and the dotted line is the 
Mecke and Dietrich resuh (Ref.[li) with Ch = 1/4. 



as the Fowler formula^^) and bending rigidity 



arei 



22.24.34. 



a.. = ldn2r'U'ir)g,{r), (86) 



sp 



32 



where U{r) is the full interaction potential. 
Since pi (z) = 0, the expression for fcgp is indepen- 
dent on the constraint used to determine pi{z). 

In Figure [21 the different models for the bend- 
ing rigidity are compared to the simulation re- 
sults of Vink et al.^ . For the evaluation of fcgp 
we have taken gi{r)=g^^ (r; -qi). 



VI. DFT: LONG-RANGED 
INTERACTIONS 

Here we examine the case that the expansion 
of a;(g, Z12) in cannot be continued to 0{q^). 
This is the case when the interaction potential 
falls of as 1/r" at large distances, with n < 6. 
In particular, we shall assume the asymptotic 
behavior of Uatt{f) to be given by 



t/att(r) = -A/r^ 



when r ^ d . 



(87) 



The analysis below only assumes that the 
asymptotic behavior of Ua.n{r) is given by the 
above expression. However, when we show ex- 
plicit results, we consider the above form for 
f/att {r) extended to the whole range 1 < r/d < 00 
(see also the Appendix). 

With the asymptotic behavior of Uatt{T) given 
by Eg. ([87]) . the expansion of uj{q, Z12) in takes 



on the form: 

tj(g, Z12) = t^o(2l2) +t^2(zi2)<7^ 

+ ^ A g4 \n{qd) -f ^^4(^12) + . • . 

The coefficient of the ln(q)-term only depends 
on the asymptotic behavior of the interaction 
potential as defined by the coefficient A, whereas 
^4(^12) depends on the interaction potential's 
full shape. With the expansion in Eq. ((88|) . a{q) 
in Eq. ()58|) can now be written in the form of 
Eq.ini)^ 

(T{q) = (T + k,q^ In(gd) -|- fco + fci q^ + 0{q'^) 
= a + ksq^ ln(g4)+0(g^) (89) 

with kg ln(£k/d) = fco + fci and 



00 00 



(90) 

ko ^ jdzi J dzi2 cd4{zi2) Po{zi)po{z2) , (91) 

— 00 —00 

00 00 

ki = Jdzi Jdzi2 u;2{zi2) Pi{zi)pq{z2) 

— 00 —00 
00 

+ ^ Jdzipiizi)p'oizi). (92) 
—00 

Next, we proceed to evaluate these expressions 
in the gradient expansion. 



A. Gradient expansion for long-ranged 
forces 

We first turn to the evaluation of kg in 
Eq. |9T|) . A straightforward gradient expansion 
of Pq{z2) is now not possible due to the fact that 
the integral Jdzi2UJ4:{zi2) is no longer finite^^-. 
The assumption of proximity to the critical 
point, however, does allow one to consider only 
the asymptotic form of uji{zi2) at large dis- 
tances. Using Eq. ([57| . one finds for ^4(212) as 
defined by Eq. 



W4(^12 



-kA 

"32" 



7e 



■In 



4^2 ) 



-O 



(93) 

One now proceeds by inserting the above expres- 
sion for L0i[zi2), together with the tanh-profile 
for po{z) in Ea. ([M)) . into the expression for k^ 
in Eq. (|9ip and carrying out the remaining inte- 
grations over zi and Z12. One finds for k^ 



ko 



32 



A{Apf 



ln(^/d)+co + 0( 



(94) 
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where 



Co = 7e - I - jdt\n{t^) 



-0.605270. 



ih(t) - tcosh(t) 



sinh''(t) 



(95) 



Next, we turn to the evaluation of ki . In the gra- 
dient expansion, the expression for fci in Eq. (j92p 
reduces to: 



ki — —2m / dz pi{z) p'q{z) 



(96) 



The further evaluation of ki requires one to 
solve the Euler-Lagrange equation in Eq. (p^ for 
Pi{z). Again, a gradient expansion of Po(z2) is 
not possible due to the fact that now the inte- 
gral J dzi2 zf2 ^^2(212) is no longer finite. Using 
the expression for the interaction potential in 
Ea. ((87)) . one finds for uj2{zi2) when |2;i2|^d 



^2(212) 



ttA 
8z^2 



^12 



(97) 



The above expression for 0^2(2:12) is used to solve 
Eq. ((62)) for pi{z) which is then inserted into the 
expression for ki in Eg. fM]) . After some algebra, 
one finally obtains for ki 



ki' 
k\' 



32 



(98) 



with 



-0.559665. 



dt \n{t) 



sinh^(t) 3sinh(t) - 3i cosh(t) 
i2 sinh^(t) 

« 1.461525... (99) 

In Figure[Tni fco , kf and kf' are shown as a func- 
tion of the reduced temperature distance to the 
critical point, t= |1 — T/Tc\. The open symbols 
are obtained from numerically solving Eas. l|49p 
and ([55]) to obtain po^z) and pi{z) and insert- 
ing the result into Eqs. ipTj) and ([5^ . The drawn 
lines are the gradient expansion approximation 
for ko in Eq.® and ki in Eq.®. 

Adding the expressions for k^ in Eq. ([94|) and 
ki in Eq. ipS)) . one obtains for 



exp(co 



\^ + 0{d) 



0.311942... C -I- 0(d), 
exp(co + cl=) e + 0{d) 
2.354329... e + C(d)- 




FIG. 10; Contributions to the bending rigidity fco 
and fci in units of k^T versus the reduced temper- 
ature distance to the critical point, t. Symbols are 
numerical results, the drawn lines are the gradient 
expansion approximation. 
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FIG. 11: The bending length Ik (in units of d) ver- 
sus the reduced temperature distance to the critical 
point, t, using the integral constraint (circles) and 
crossing constraint (squares) . Open symbols are nu- 
merical results, the drawn lines are the gradient ex- 
pansion approximation. The dashed line is the cor- 
relation length ^ and the dotted line is the Mecke 
and Dietrich result (Ref. [l^ with Ch = 1/4. 



In Figure [TTl the bending length is shown as a 
function of t. The open symbols are the numer- 
ical results. The drawn lines are the gradient 
expansion approximations in Eq. ljlOOp . 

Other than that the gradient expansion seems 
not to be as accurate in reproducing numerical 
results, the results in Figures [TO] and [11] are in 
line with the earlier results obtained for short- 
ranged forces. The leading order correction to 
the surface tension is negative when ql^ < 1, 
and the effect is more pronounced for the cross- 
ing constraint than the integral constraint. Al- 
though the goal of Mecke and Dietrich in ref. [l^ 
is to include higher order terms, terms beyond 
in the expansion of cr(g), it is also interest- 
ing to compare with the Mecke and Dietrich ap- 
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proach for the terms obtained to order ln('?) 
and . Using the Mecke and Dietrich expression 
for pi{z) in Eq. (|50p . one has to leading order in 
the gradient expansion: 

CTMD(g) =(T + ks q^ \\i{qd) + koq'^ + kuB q^ + ■■ 
= a + ksq^Hq£f^) + 0{q^), (101) 

where fcg and ko are given by the expressions in 
Eg. ((92)) and Ea. ((94)) . and where kuD is given by 
the previously derived expression in Eq. (|85p for 
short-ranged forces. In Figure [11] we show, as 
the dotted line, the result for the bending length 

VII. DISCUSSION 

In the first part of this article, we have demon- 
strated that the full spectrum of surface fluc- 
tuations obtained in Monte Carlo simulations^ 
of the colloid-polymer interface, can very accu- 
rately be described by the following expression: 

S(q) = ^-'l^+MLS,{q). (102) 
(T q- o"^ 

The three terms in this expression work in three 
different g-regimes: 

1. Classical capillary wave regime, gd<Cl 

2. Extended capillary wave regime, qd<l 

3. Bulk-like fluctuations regime, qd>l 

Two adjustable parameters are present: A/l that 
weighs the bulk-like fluctuations compared to 
the capillary wave fluctuations and the bend- 
ing rigidity k as defined by the leading or- 
der correction in an expansion of <j{q) in q^, 
(T{q) = a + kq^ + . . .. We found that a fit to 
the simulation results yields fc<0, i.e. the lead- 
ing order curvature correction tends to lower the 
surface tension (j{q). This effect is termed cap- 
illary enhancement^; capillary waves are "more 
violent", less restricted by surface tension, at 
smaller wavelengths. 

One could worry whether a negative bending 
rigidity is consistent with having a stable inter- 
face. Tarazona et al.'^ indicate that a decrease of 
(j{q) = a+k q^+- ■ ■ ultimately leads to a destabil- 
isation of the interface at large q. It is therefore 
important to realise that the extension of the 
capillary wave model, through the inclusion of a 
bending rigidity, is valid only for low q. Higher 
order terms in the expansion in q are not sys- 
tematically included. In this article we propose 
to describe S{q) for large q {qd> 1) in terms of 
molecular, bulk-like fluctuations through Sb{q)- 
It is shown that the full S{q), which is then 
a combination of the extended capillary wave 



model at low q and bulk-like fluctuations at large 
q, remains well-behaved ensuring the stability of 
the interface. For systems with a low (or even 
zero) surface tension, the bending rigidity is the 
dominant contribution near q ~ and one nec- 
essarily requires a positive value for k^^ , but for 
the simple, (quasi) one-component system con- 
sidered here this is not an issue. 

A most important and generally underappre- 
ciated point that we like to emphasize is that the 
location of the interface cannot be deflned un- 
ambiguously. A certain procedure must always 
be formulated to determine the height function 
/i(fjj). We have shown that different choices 
for the location of the interface, which are all 
equally legitimate as long as they lead to a lo- 
cation of the dividing surface that is 'sensibly 
coincident' with the interfacial region^^, lead to 
different results for the bending correction to 
the capillary wave model. Naturally, all experi- 
mentally measurable quantities cannot depend 
on the chosen location of the interface, mak- 
ing it necessary to formulate precisely the quan- 
tity that is determined in experiments or sim- 
ulations. It was shown that for the simulation 
results, the value of the bending rigidity in the 
above expression for S{q), corresponds to the 
height function being defined according to the 
integral constraint, k = k'^'^. 

For the determination of fc, it is necessary to 
take the contribution from bulk-like fluctuations 
into account since they also contribute as a con- 
stant, A/lS'6(0), in the capillary wave regime 
{qd^ 1). This observation is consistent with the 
interpretation of light scattering results by Dail- 
lant and coworkers^. To determine a{q), they 
subtract from S{q) a contribution proportional 
to the penetration depth (cx A/l) times the liq- 
uid compressibility (cx Sb{Oj)- Even though the 
light scattering results by Daillanti are obtained 
for real fluids, for which the interaction poten- 
tial is not necessarily short-ranged, one expects 
that a description in terms of the above men- 
tioned three regimes is again useful. A further 
comparison with the light scattering results is, 
however, necessary. 

In the second part of this article, a molecular 
theory to describe the inclusion of the bending 
rigidity correction to the capillary wave model is 
presented. An essential feature of the theory is 
the 'Ansatz' made in Eg. ([55)) regarding the ther- 
modynamic conditions used to vary the interfa- 
cial curvature. It improves on earlier choices 
made in the sense that the bulk densities are 
equal to those at coexistence and the density 
profile is a continuous function^ii^i^ii^S. The 
theory predicts that the scaling behavior of the 
bending rigidity equals that of the surface ten- 
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sion near the critical point 

k (X (T d'^ (X , 



(103) 



where /i« 1.26 is the usual surface tension criti- 
cal exponent (in mean- field /i = 3/2)''. This new 
scaling prediction differs fundamentally from the 
scaling of the bending rigidity in the 'equilib- 
rium approach', fcoqoccr^^. In this approach the 
bending rigidity is determined from considering 
the equilibrium free energy of spherically and 
cylindrically shaped liquid droplets, with their 
radii varied by changing the value of the sys- 
tem's chemical potentia P°'^^i'^^ . 

The negative sign and scaling behavior of the 
bending rigidity obtained from the molecular 
theory are in accord with Monte Carlo simu- 
lations. However, the magnitude of k from the 
molecular theory, y/—k/a ~ 0.13 d, is signifi- 
cantly below the value obtained in the simula- 
tions, y/—k/a w 0.47 d (see also Figure [8]), but 
we believe this to be due to simplifications made 
in the theory rather than a true discrepancy. 
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APPENDIX A: DEPLETION 
INTERACTION POTENTIAL 

The phase-separated colloid-polymer system 
is effectively treated as a one-component system 
considering the colloids only. The colloid-colloid 
interaction is then given by a hard sphere repul- 
sion (diameter d) with an attractive depletion 
interaction^^ induced by the presence of poly- 
mers (radius i?g): 



Udepir) 



-ksTrip 



2(e-l)3 



2e-^ -3e 



'r\3 
.d) 
(Al) 

with d<r <ed and the size ratio parameter e is 
defined as: 

2R, 
~d 

Using this form for Uattir), the coefficients a, 
m, B are readily calculated to yield 

a = kBTd^r]p^{2 + 6e + 3e^ +e^), 



1 



(A2) 



m = kBTd^T]p^{5 + 15e + lQe^ + 6e^ 



h3e4 + £5) 



(A3) 



B = kBT d^ Tjp -^{28 + 8Ae + 63e^ 

-h45 + 30 £4 + 18 ^5 _^ 9 £0 + 3 e'^) . 

One may also determine the functions 1^0(212), 
W2{zi2) and W4(zi2). When |zi2|<fi one has: 
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When d < I Z12 1 < £(i one has: 
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APPENDIX B: LONDON-DISPERSION 
FORCES 

The following explicit form for Uatt (f ) is con- 
sidered: 

C/att(r) = I when r>d, ^^^^ 

I when r < d . 

Using this form for ?7att('"), the coefficients a and 
m are readily calculated to yield 

2'n-A , ttA 

a — and m — — - . (B2) 



3d3 



3d 



18 



With this form for the interaction potential one 
may expand ui{q, 212). When \zi2\ > d one has: 



TtA 4 
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When |2;i2| <rf one has: 
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